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In [1] the authors studied the spectrum of the Chebyshev collocation 
operator for the heat equation, proving that it shares certain 
properties with the spectrum of the corresponding differential 

operator. Here we present similar results for spectral discretization 
of hyperbolic systems. 

Assume that the differential system, under given boundary conditions 
has only solutions which decay in time to zero. To test whether the 
numerical solutions behave in the same way, it is sufficient to check 
the corresponding scalar problem with zero inflow: 

U t = u x ln l x l < 1, t > 0 

(1) ) 

u = 0 at x = 1 

If all the numerical solutions of (1) decay, so will the solutions of 
the general system, provided that the spectral approximation contains an 
even or odd forcing polynomial (to be subsequently defined). 

In order to prove this statement we require the following two 
lemmas : 


Le ““ '• isi. P (x) be a real Hut wltz polynomial (l.e., having only 
roots with negative real parte), end let e be a complex -...her, less 
than one in magnitude. Then the polynomial 


<i(x) = p(x) + ep(-x) 


is also Hurwitz. 



Proof. Let us first recapitulate a few facts about positive 
pairs. Two polynomials g(x) and h(x) are said to form a positive 

pair if: 

a. they have positive coefficients and real, negative and distinct 
roots 

b. their degrees are equal, or degree(h) = degree(g) -1 

c. their roots interlace: 

0 > uj > Vj > U£ > V£ etc. , 

where u ± are the roots of g and v^ the roots of h. The basic 
theorem states that the polynomial f(x) is Hurwitz, if and only if, 

f(x) = g(x 2 ) + xh(x 2 ) 
with g and h forming a positive pair. 

Without loss of generality, let p have real positive coefficients; 
write p using a positive pair as: 

p(x) = g(x 2 ) + xh(x 2 ). 

2 

The polynomial q is Hurwitz, if and only if, Q(x) = |q(x)| is. 
For Q we have the expression: 

Q(x) = (p(x) + ep(-x)) (p(x) + e p(-x)) 

2 2 2 2 

= (1 + 2 Re e + |e|^)g^(x^) + (1 - 2 Re e + |e| )x h (x ) 

+ x(2(l - |e| 2 )g(x 2 )h(x 2 )) 
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and all its coefficients are positive. Thus we must prove that: 

G(x) = (1 + 2Re e + |e| 2 )g 2 (x) + (1 - 2Re e + |e| 2 )xh 2 (x) 

and 

H(x) = 2(1 - | e| 2 )g(x)h(x) 

form a positive pair. After simplifying positive factors, it is enough 
to show that 

G (x) = A 2 g 2 (x) + xh 2 (x) 

and 

H + (x) = g(x)h(x) 

form a positive pair, where A 2 = (1 + 2Re e + |e| 2 )/(l - 2Re e + |e| 2 ). 
As g and h form a positive pair and A is positive, the polynomial 
Ag(x 2 ) + xh(x 2 ) is Hurwitz, and so is its square: 

( Ag(x 2 ) + xh(x 2 )) 2 = A 2 g 2 (x 2 ) + x 2 h 2 (x) + x[2Ag( x 2 )h(x 2 )). 

By examining the last formula we recognize G + and H + , which form a 
positive pair, as required. 
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T-Pimna 2. Let <|>(x), B(x) be polynomials of degree M. Then the 
solution of 

( u = ± u + a(t)B(x) in t > 0, |x| < » 

t t 

( 2 ) < 

I u(x,t=0) = 4>(x) 

is 

t 

(3) u(x, t) = <Kx ± t) + / o(9)B(x ± 9)d0. 

0 

Proof. Take the Laplace transform 


Then: 


f(t) * f(x) = / e Tt f(t)dt. 

0 


TU - <|> = ± -|^ + aB 


Since <{• and B are polynomials in x one may invert formally the 
£ 

operators (x ± to obtain: 


u(x,x) = l 
k=0 


(±l) k B (k) (x) - 


k+1 


a(x) + l 
k=0 


(±l) k 4> (k) (x) 
k+1 


= l a(x) 


B (k) (x) 

k! 


(± 1 ) 


k k! 
k+1 


k=0 


<t» (k) (x) 
k! 


(±D 


k kl 
k+1 


As -yiy is the Laplace transform of t k , we recognize the transform of 
x k+ 

the Taylor series for B(x ± t), 4» Cx ± t). We can write: 
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u(x,t) = a B(x ± t) + «))(x ± t) 

and finally obtain (3) by interpreting the product above as the 
transform of a convolution. Remark that the solution u(x,t) is itself 
a polynomial of degree M in x, with time dependent coefficients. 
Consider now the first order hyperbolic system: 


(5) 



u 


t 


u 

X 



u = Rv 


v = Lu 


in | x | < 1, t > 0 


at x = 1 


at x = -1 


Here u, v are column vectors of dimensions m and n 
respectively. These variables are coupled only in the boundary 
conditions via the two constant matrices L and R, of dimensions 
m x n and n x m. 

We shall study this system and its spectral discretizations to prove 
that if all the solutions of (5) tend to zero for large times, so do the 
numerical solutions. 

First, remark that there is decay in time, if and only if, both 
matrices LR and RL are contractions, i.e., 
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BLRx# < Ixll , lIRLxll < llxll 


for all x * 0, in some suitable norm. This is obvious if L and R 
are interpreted as reflection coefficients, since u and v are 

transported without change on characteristics between x = 1 and 
x = -1. 

Next we introduce the spectral approximation for system (5). Define 
the Chebyshev points 

x. = cos(nj/N) 0 < j < N 

J _ - 


and collocate at the interior points x i » x 2* * * » x n- 1* This means that 

the numerical solutions u N ,v N are polynomials in x of degree < N, 
subject to the boundary conditions (5c-d), which also satisfy the 
differential equations (5a-b) at x = Xj , 0 < j < N. Since the 
polynomial T^(x) vanishes at exactly these points and is of degree 
N - 1, we find that the spectral solutions satisfy: 


( 6 ) 





9u N 

— + “ (t) V x) 

8v m 

- 7T + 6(t)T N (x) 


“n “ Ev n at x = 1 




= Lu 


N 


at x = -1 
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(we assume that the time evolution is not discretized). Using Lemma 2 
we can solve (6a-b) in the form: 


<j>(l+t) + a T'(l+t) = R(iKl-t) + 3 T'(l-t)) 


♦ ( 


-1-t) + $ T'(-l-t) = L(<K-l+t) + a T'(-l+t)) 


or 


i ra+t) -r r(i-t) 

m N N 


-l r(-i+t) i r(-i-t) 

N n N 




r_- 

a 


a 


= A 



def 






RK l-t) - <Ki+t) 


L<f> ( — 1+t ) + <(>(-1— t) 


def 


(I m and I n are identity matrices of the corresponding dimensions.) 
Since A and r are rational functions of x (as transforms of 
polynomials) it is clear that a and F will be rational. This means 
that the original functions ct,B must be of the form: 


(7) 


l 


V 

e p k (t) 


for some polynomials p k ; moreover the exponents are precisely the 
poles of a and 3, These poles may be values where det | A [ =0, or 
poles of the right-hand-side. However, we have the expansion: 


A = (A o + A i T + 

T 
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r = -M U 0 + r l 

T 


with M < N, since the degree of the coefficients is exactly N - 1, 
while the degree of r is < N. Thus x = 0, which is the only pole of 
the right— hand— side cannot be a pole of a or B and we only have to 
discuss the zeros of det | A | « Finally, since a function of the form 
(7) decays, if and only if, ReX fc < 0 for all k, we are left to prove 
that the determinant of A vanishes only in the open left half of the 
complex plane. 

Introduce the polynomials of degree N — 1 p(x) and q( T )> defined 


Tju+Ey = — P (T) 

T 


T N (1 " t) = “N q(T) * 
T 


It is easily checked that q(x) = p(-x)(-l) 5 moreover T' is even or 

odd, satisfying 


T-(l) = (-1) N_1 T'(-l) . 


Using these symmetry properties, det |a| =0 implies that 


I p(t) -Rp(-x) 


-Lp(-x) I p(x) 
n 
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Then, the matrix identity 


I m P(T) -Rp(-T)\/l m p(x) Rp ( —X ) ^ 


I P 2 (t) 
m 


-Lp(-x) I_p(x)/ \ 0 I n p(x)/ \-L P (t) p (-t) I n p 2 (x) - LRp 2 (— x) y 


shows that the solutions of (9) may be found among the solutions of 
p(x) = 0 and the solutions of 

det |p 2 (x)I n - p 2 (-x)LR| = 0. 

This last equation simplifies to: 


( 10 ) 


p(x) ± Xp(-x) = 0 


where X is any eigenvalue of LR. We know that |x| < 1, since LR 

is a contraction. Moreover, in [2] it was shown that the polynomial 

p(x) is Hurwitz. Using Lemma 1, we conclude that the solutions x of 
(10) have negative real parts, and a and g tend to zero at large 
time, decreasing exponentially. This is sufficient to make all 

solutions of (6) vanish as t - > «. Indeed, let u and v solve (5) 
with polynomials of degree < N as initial values. Then the 

differences u " u N > v - v N satisfy the equation (2) with <(. « 0, and 

the explicit formula (3) shows they decay; since u and v also decay, 
by assumption, our assertion is proved. 


T 
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The proof we presented hinges on two facts: 

a. p(x ) is Hurwitz - exactly the required property for the 
scalar equation (1). 

b. formula (8) - a symmetry property of T' . 

Many other spectral approximations are possible — using different 
sets of orthogonal polynomials, or Galerkin and x— methods instead of 
collocation. All of them will replace (5) by an equation similar to 
(6), the only difference being that T' must be changed to another 
polynomial (e.g., % for the x- method). This polynomial is known 
explicitly for any given spectral method and we shall call it the 
forcing polynomial, since it appears as an inhomogenous term in (6). If 
the forcing polynomial has the correct symmetry and the corresponding 
scalar problem has only decaying solutions, the method of proof readily 
extends to the general spectral method. We should also remark at this 
point that the most usual Chebyshev discretization for u t = u x , 
collocating at the interior points and at the outflow, has the forcing 
polynomial (x+l)T' (x) , which is neither odd nor even, and thus is not 
covered by our theory. 

In conclusion, we summarize our results in the following 

Theorem. Let the solutions of (5) decay to zero in time. Consider 
a spectral approximation for (5) possessing an even or odd forcing 
polynomial, such that the numerical solutions obtained by this method 
applied to (1) tend to zero as time increases. Then the numerical 
solution corresponding to the system (5) also tends to zero as time 


increases. 
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